using System;

namespace BOL.Maths.LinearAlgebra
{
	/// <summary>
    /// LU Decomposition.
    /// This method decomposes the matrix A according to L * U = A
    /// where L is a lower triangular matrix and U is an upper triangular matrix
	/// </summary>
	public sealed class LUDecomposition : DecompositionBase
	{
		#region Private variables

        private double[][] _lu; // matrix for internal storage of decomposition.
        private int _pivotSign;
        private int[] _pivot;
		
		#endregion

		#region Public properties

        /// <summary>Returns the determinant.</summary>
        public double Determinant
        {
            get
            {
                if (_m != _n)
                    throw new ArgumentException("Matrix must be square.");

                double det = (double)_pivotSign;
                for (int i = 0; i < _m; i++)
                    det *= _lu[i][i];

                return det;
            }
        }

		/// <summary>Returns whether the matrix is singular or not.</summary>
		public bool IsSingular
		{
			get
			{
                for (int i = 0; i < _m; i++)
                    if (_lu[i][i] == 0)
						return true;

				return false;
			}
		}

		/// <summary>Returns a lower triangular matrix.</summary>
		public double[][] L
		{
			get
			{
                double[][] lower = new double[_m][];
                int i, j;

                for (i = 0; i < _m; i++)
                {
                    lower[i] = new double[_n];
                    for (j = 0; j < _n; j++)
                    {
                        if (i > j)
                            lower[i][j] = _lu[i][j];
                        else if (i == j)
                            lower[i][j] = 1.0;
                        else
                            lower[i][j] = 0.0;
                    }
                }

                return lower;
			}
		}

		/// <summary>Returns a upper triangular matrix.</summary>
		public double[][] U
		{
			get
			{
                double[][] upper = new double[_n][];
                int i, j;

                for (i = 0; i < _n; i++)
                {
                    upper[i] = new double[_n];
                    for (j = 0; j < _n; j++)
                    {
                        if (i <= j)
                            upper[i][j] = _lu[i][j];
                        else
                            upper[i][j] = 0.0;
                    }
                }

                return upper;
			}
		}

		/// <summary>Returns _pivot permutation vector</summary>
        public double[] PivotPermutationVector
        {
            get
            {
                double[] temp = new double[_m];

                for (int i = 0; i < _m; i++)
                    temp[i] = (double)_pivot[i];

                return temp;
            }
        }

		#endregion
		
		#region Constructor
		
		/// <summary>LU Decomposition</summary>
		/// <param name="matrix">a rectangular matrix</param>
		public LUDecomposition(double[][] matrix) : base(matrix)
		{
            _lu = new double[_m][];
            _pivot = new int[_m];
            for (int i = 0; i < _m; i++)
            {
                _lu[i] = new double[_m];
                Array.Copy(matrix[i], _lu[i], _n);
                _pivot[i] = i;
            }
			_pivotSign = 1;

            InitializeLU();
		}

        private void InitializeLU()
        {
            int i, j, k;
            double[] row;
            double[] column = new double[_m];
            int kmax;
            double sum;

            for (j = 0; j < _n; j++)
            {
                // Make a copy of the j-th column to localize references.
                for (i = 0; i < _m; i++)
                    column[i] = _lu[i][j];

                // Apply previous transformations.
                for (i = 0; i < _m; i++)
                {
                    row = _lu[i];

                    // Most of the time is spent in the following dot product.
                    kmax = Math.Min(i, j);
                    sum = 0.0;

                    for (k = 0; k < kmax; k++)
                        sum += row[k] * column[k];

                    row[j] = column[i] -= sum;
                }

                // Find _pivot and exchange if necessary.
                int p = j;
                double t;
                int k2;

                for (i = j + 1; i < _m; i++)
                    if (Math.Abs(column[i]) > Math.Abs(column[p]))
                        p = i;

                if (p != j)
                {
                    for (k = 0; k < _n; k++)
                    {
                        t = _lu[p][k];
                        _lu[p][k] = _lu[j][k];
                        _lu[j][k] = t;
                    }

                    k2 = _pivot[p]; _pivot[p] = _pivot[j]; _pivot[j] = k2;
                    _pivotSign = -_pivotSign;
                }

                // Compute multipliers.
                if (j < _m & _lu[j][j] != 0.0)
                    for (i = j + 1; i < _m; i++)
                        _lu[i][j] /= _lu[j][j];
            }
        }

		#endregion

		#region Public Methods

        /// <summary>Solve A * X = B</summary>
        /// <param name="b">A Matrix with as many rows as A and any number of columns.</param>
        /// <returns>X so that L * U * X = B(piv,:)</returns>
		public double[][] Solve(double[][] b)
		{
            if (b.Length != _m)
				throw new ArgumentException("Matrix row dimensions must agree.");

			if (IsSingular)
				throw new SystemException("Matrix is singular.");

            int i, j, k; // counts for row, column, 

			// Copy right hand side with _pivoting
			int nb = b[0].Length;
			double[][] x = MatrixDouble.SubMatrix(b, _pivot, 0, nb - 1);
			
			// Solve L * Y = B(_pivot,:)
            for (k = 0; k < _n; k++)
                for (i = k + 1; i < _n; i++)
                    for (j = 0; j < nb; j++)
                        x[i][j] -= x[k][j] * _lu[i][k];

			// Solve U * X = Y;
            for (k = _n - 1; k >= 0; k--)
			{
                for (j = 0; j < nb; j++)
                    x[k][j] /= _lu[k][k];

				for (i = 0; i < k; i++)
                    for (j = 0; j < nb; j++)
                        x[i][j] -= x[k][j] * _lu[i][k];
			}

			return x;
		}

		#endregion
	}
}
